email: mperez.mendez@izt.uam.mx
url: https://github.com/MarcoAPerezM
url_lab: https://sites.google.com/view/deserlab/
El análisis de series de tiempo se configura como una herramienta fundamental entre el conjunto de herramientas con las que cuenta el economista moderno. El mercado laboral del economista demanda el uso de instrumentos analíticos que le permita realizar un diagnóstico preciso sobre los diferentes fenómenos económicos, políticos, sociales y financieros. En los tiempos modernos, las herramientas se enfocan, cada vez más, en un uso intensivo de datos acompañado de capacidades analíticas derivadas de la aplicación de la teoría económica. Entre estas se encuentra el modelado matemático de fenómenos y su correspondiente contrastación empírica. El análisis de series de tiempo está, entre las herramientas econométricas, ubicado como un conjunto de instrumentos que el mercado laboral demanda y valora en demasía. En este libro se encuentra el desarrollo de las series de tiempo univariadas. En el capítulo 1 se desarrollan las características de las series de tiempo, las herramientas básicas para su análisis y se esgrime la necesidad de desarrollar modelos sofisticados, pero sencillos, que permitan entender el pasado de los fenómenos y predecir su futuro. En el capítulo 2 se desarrolla el Modelo de Componentes Subyacentes: tendencia, ciclo y estacionalidad se modelan desde la perspectiva clásica. En el capítulo 3 se desarrollan filtros y suavizamientos. Medias móviles simples, exponenciales, etc. muestran como se comportan los fenómenos en el largo plazo. Se muestra como el filtro Holt-Winters permite reproducir con un alto grado de precisión los fenómenos con tendencia y estacionalidad y, además, permite obtener pronósticos altamente eficientes. En el capítulo 4 se desarrolla la metodología Box-Jenkins y se construyen los modelos tipo SARIMA.
Como es bien sabido, todo Análisis ecónomico de series de tiempo requiere de ciertos prerequisitos imprescindibles para poder desarollar el análisis en cuestión. En particular, este libro requiere el conocimiento y dominio de diferentes tópicos que le son imprescindibles al economista moderno.
En primera instancia, se encuentra el entendimiento de herramientas estadísticas y de probabilidad que se requieren, tanto para el desarrollo de instrumentos como de indicadores que permitan reproducir los patrones de comportamiento subyacente de los fenómenos económicos, al mismo tiempo sirven para evaluar el grado de exactitud de sus reproducciones.
También se requiere el entedimiento, en un nivel básico, de teoría econométrica: estimación puntual, por intervalos, por cuantiles, etc. Técnicas de estimación de parámetros y modelado matemático.
Otro de los requisitos fundamentales es saber identificar con precisión las fuentes de datos donde se encuentran las series de tiempo relevantes para el análisis económico.
Cuando el profesional cuenta con datos, tiene dos opciones: 1) la primera es cargar directamente en formato csv (Variables Separadas por Comas) y la segunda es cargarglos en formato xlsx (Excel).
Para abrir datos en csv se debe emplear la función read.csv:
Para abrir datos en Excel se debe emplear la función read_excel() perteneciente a la librería “readxl”, así como indicar la hoja en la que se encuentran los datos:
# install.packages("readxl")
#library(readxl)
#datos2 <-read_excel("/dirección/libro.xlsx", sheet="hoja")
#datos2Una manera mucho mas eficiente y moderna de obtener datos es descargarlos directamente desde repositorios oficiales en donde se encuentran los datos. Por ejemplo, la libreria Quantmod (https://www.quantmod.com. Quantitative Financial Modelling and Trading Framework for R) brinda acceso a las bases de datos de Google, Yahoo, oanda y la FRED (Federal Reserve of Economic Data).
Por ejemplo, por medio de la librería getFX() es posible descargar el tipo de cambio del peso mexicano contra el dolar estadounidense desde oanda:
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.3.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## [1] "USD/MXN"
o el rublo ruso contra el dolar estadounidense:
## [1] "USD/RUB"
o el precio del oro en relación con el dolar estadounidense:
## [1] "XAU/USD"
Con la función getSymbols() se puede desarcargar el Bitcoin desde la FRED, sólo basta con identificar la clave de la variable de interés, en este caso “CBBTCUSD”:
## [1] "CBBTCUSD"
En este caso, el interés principal se ecuentra en analizar series de tiempo mexicanas. Para ello, es necesario conectarse a la API de INEGI (Instituto Nacional de Geografía e Informática https://www.inegi.org.mx/servicios/api_indicadores.html). Una API (Interfaz de Programación de Aplicaciones) es, en términos generales, un conjunto de protocolos o reglas que permiten la interacción entre diferentes aplicaciones de software con el objetivo de que se comuniquen entre sí para intercambiar, principalmente, datos, características y funcionalidades. Para tener acceso a la API de INEGI es necesario solicitar una llave de acceso conocida como token.
Existe una librería que permite la conexión entre R y la API de INEGI.
Si es de nuestro interés analizar el comportamiento del Poducto Interno Bruto Trimestral, es necesario identificar su ubicación en el Banco de Información Económica. La ubicación exacta es Indicadores económicos de coyuntura > Producto interno bruto trimestral, base 2018 > Series Originales > Valores a precios de 2018. EL identificador del PIB trimestral es “735879”. La función para realizar esta conexión es inegi_series(). Es necesario indicar el identificador de la variable, el token y la base de datos, la cual puede ser el Banco de Indicadores o el Banco de Información Económica.
Como puede observarse, por medio de la función plot(), la descarga se realiza con el orden inverso, la tabla se encuentra ordenada de los valores más actuales a los más antiguos.
Es necesario verificar la primer fecha de publicación de la variable y la frecuencia con la que se repite dentro de un año.
## date date_shortcut values notes
## 174 1981-04-01 Q2 11564025 NA
## 175 1981-01-01 Q1 11345848 NA
## 176 1980-10-01 Q4 10927666 NA
## 177 1980-07-01 Q3 10392733 NA
## 178 1980-04-01 Q2 10342350 NA
## 179 1980-01-01 Q1 10401368 NA
Como puede observarse, la fecha de inicio es el primer trimestre de 1980. Es necesario invertir el orden, para ello diseñamos una función que cambie el orden de acuerdo a la fecha y mantenga solo la columna de valores. A esta función personalizada la hemos llamado ts_invert().
Al aplicar esta función personalzida sobre la variable pib tenemos el orden adecuado
Ahora es necesario transformar la variable en una serie de tiempo. Para ello usamos la función ts() la cual permite identificar la fecha de inico y la frecuencia. En este caso, la variable pib se transformará en una serie de tiempo que inicia en el primer trimestre de 1980 y tiene frecuencia 4, ya que es trimestral.
Al momento de realizasr un análisis de series de tiempo siempre se va a cometer algún tipo de error. El más frecuente es el error de predicción o de estimación, este se define como la diferencia entre el valor observado y el valor estimado:
\[\hat{e}_t = y_{t} − \hat{y}_t \]
Estos errores son causados por un gran número de fuentes, por ejemplo, hay errores de medición en los procesos de recolección de las variables, hay errores de muestreo, hay errores derivados de la aleatoriedad del comportamiento humano, etc. Estos errores pueden ser grandes y sistemáticos y esto depende, en gran medida de:
a) Identificación errónea de los patrones subyacentes y relaciones entre variables. Podemos identificar, de manera equivocada, un patrón equivodado. También podemos construir un modelo estadístico basado en pocos datos y suponer que el patrón es estable a lo largo de mucho tiempo.
b) Patrones inexactos. En general, en las ciencias sociales y, en particular, en la Economía los patronos que podemos identificar son, por naturaleza, inexactos e imprecisos. Esto se debe a que no es una ciencia determinista como la física o la quimica. Las relaciones económicas no son deterministas, son aleatorias y a pesar de que podemos identificar patrones de comportamiento estables en el tiempo, siempre habrá incertidumbre y múltiples factores que generan variaciones en los resultados. Por ejemplo, podemos construir un modelo que capture, relativamente bien, el ingreso de los hogares en México. Sin embargo, ¿cuántas variables se encuentran involucradas en el la determinación del ingreso de un hogar particular? Supongamos que el ingreso promedio de un hogar es de $12,000. ¿Cuántos hogares ganan esa cantidad? De dichos hogares, ¿cuántos tiene jefatura femenina? ¿Cuántos jefes de hogar tienen educación universitaria? etc. etc. Hay un gran número de variables involucradas en la determinación del ingreso, todas esas variables generan una fuerte variabilidad e impiden que el fenómeno sea determinista, por ello es imposible identificar con exactitud el ingreso de los hogares con base en un modelo estadístico o econométrico. Esta variación induce un error de estimación.
c) Cambios en los patrones. Siguiendo el ejemplo del ingreso de los hogares, un error común es suponer que los patrones que se identificaron en un momento en el tiempo serán los mismos durante un periodo de largo plazo. Esto no es necesariamente cierto, hay cambios en la estructura de los hogares, en la estructura productivia, en los niveles de educación y productividad, etc. Estos cambios inducen cambios en los patrones subyacentes de los fenómenos y por lo tanto en los parámetros de los modelos, a este dilema normalmente se le conoce como cambio estructural o paramétrico.
Por ello es necesario desarrollar herramientas que nos permitan identificar el grado de precisión de los modelos construidos. Para evaluar la precisión vamos a suponer que tenemos un fenómeno \(y_t\) observado.
Vamos a suponer que tenemos dos modelos alternativos para reproducir y, eventualmente, predecir el fenómeno \(y_t\).
yp1 <- c(23,25,26,29,31,34,41,40,38,44)
yp2 <- c(24,26,28,29,33,34,45,38,43,45)
y <- ts(y)
yp1 <- ts(yp1)
yp2 <- ts(yp2)
ts.plot(y,yp1,yp2, col=1:3)\[ EMA =\frac{1}{T} \Sigma_{t=1}^T | e_t |\] El EMA es una medida de precisión general que brinda una idea del grado de dispersión y cuenta con la característica de que le brinda la misma importancia a todos los errores.
## [1] 3.666667
Construimos funciones personalizadas para cada indicador.
## [1] 3.666667
## [1] 3.333333
Es una medida de precisión que indica que tan dis perso está el modelo de los datos reales pero cuenta con la característica de brindarle mayor peso a los errores más grandes
\[ ECM= \frac{1}{T}\Sigma_{t=1}^T(e_t)^2 \]
## [1] 21.88889
## [1] 20.88889
Es la medida relativa del EMA
\[ EMAP= \frac{1}{T}\Sigma_{t=1}^T\frac{|e_t|}{y_t}\times100 \]
## [1] 9.856951
## [1] 8.889399
\[ ECMP= \sqrt{\frac{1}{T}\sum_{t=1}^T \left(\frac{e_t}{y_t} \right)^2}\times 100 \]
## [1] 22.36365
## [1] 9.784647
\[U = \sqrt{\frac{\frac{1}{T}\sum_{t=1}^Te_t^2}{\frac{1}{T}\sum_{t=1}^Ty_t^2}} \]
## [1] 0.08012821
## [1] 0.03846154
ts_precision <- function(x, y,
type = c("ecm", "ecmp", "ema", "emap",
"UTheil"),
na.rm = TRUE)
{
switch(match.arg(type),
ecm = ecm(x),
ema = ema(x),
emap = emap(x, y),
ecmp = ecmp(x, y),
UTheil=UTheil(x,y))
}
ema <- function(x){
l <- length(x)
s <- sum(abs(x))
s/l
}
ecm <- function(x){
l <- length(x)
s <- sum(x^2)
s/l
}
emap <- function(x,y){
l <- length(x)
s <- sum(abs(x)/y)*100
s/l
}
ecmp <- function(x,y){
l <- length(x)
s <- (sum(x/y)^2)
ss <- sqrt(s/l)*100
ss
}
UTheil <- function(x,y){
s <- (sum(x)^2)/(sum(y)^2)
ss <- sqrt(s)
ss
}
ts_precision(error1,y,type="ecm")## [1] 21.88889
Construimos una función que además de calcular todos los indicadores de precisión, los pone en una tabla comparativa.
ts_precision_all <- function(x,y){
ema <- ts_precision(x,y,type="ema")
emap<- ts_precision(x,y,type="emap")
ecm <- ts_precision(x,y,type="ecm")
ecmp <- ts_precision(x,y,type="ecmp")
u <- ts_precision(x,y,type="UTheil")
tabla <- as.matrix(data.frame(ema, emap, ecm, ecmp, u))
colnames(tabla)<- c("EMA", "EMAP", "ECM", "ECMP", "UTheil")
knitr::kable(
tabla, booktabs = TRUE,
caption = "Indicadores de Precisión"
)
}| EMA | EMAP | ECM | ECMP | UTheil |
|---|---|---|---|---|
| 3.666667 | 9.856951 | 21.88889 | 22.36365 | 0.0801282 |
| EMA | EMAP | ECM | ECMP | UTheil |
|---|---|---|---|---|
| 3.333333 | 8.889399 | 20.88889 | 9.784647 | 0.0384615 |
Los modelos de series de tiempo de Patrones de Comportamiento Subyacente pretenden recoger las regularidades de la serie capturadas en el comportamiento de una variable a lo largo del tiempo. Descansab en la idea de que los componentes subyacentes se pueden reproducir y nos permiten analizar los datos observados. La metodología es muy básica: se observa un determinado fenóomeno, se analizan sus regularidades o patrones y se construyen modelos basados en dichas regularidades observadas. Descansan en la idea de que toda serie de tiempo se puede descomponer en diferentes elementos subyacentes que forman parte de la serie, de manera indirecta, y que pueden explicar su evolución, a lo largo del tiempo. Los Modelos del Componentes Subyacentes (MCS) no pretenden representar el proceso generador de datos sino, mas bien, describir, de manera muy general, las características fundamentales de las series.
La construcción de que una serie de tiempo se puede construir por elementos superpuestos no es nueva, ya en la edad media se preguntaban por que había fenómenos qe se repetían de manera regular en el tiempo. El punto de partida para la construcción de este tipo de modelos es, por lo tanto, saber identificar la existencia de dichos patrones o regularidades. Este enfoque pretende clasificar los tipos de movimientos que caracterizan una serie de tiempo como tendencia, estacionalidad, ciclo e irregular. Estos cuatro componentes son considerados como los patrones de comportamiento subyacente.
El objetivo del análisis de series de tiempo es dentificar los patrones subyacentes expresados en un modelo uniecuacional.
\[Y_t=f(T_t, C_t, E_t, I_t) \] El componente irregular \(I_t\) representa todos los movimientos no sistmáticos de la serie a lo largo del tiempo. Son las perturbaciones aleatorias que se presentan en la vida cotidiana. Son impredecibles y en suma su efecto se anula entre ellas. La estacionalidad es un patrón que se reproduce al interior de un año y refleja los impulsos del fenómeno asociados con efectos observados en diversas estaciones del año. El ciclo es un comportamiento de mediano plazo que, normalmente son obserbables en series con temporalidades superiores a un año.Por último, la tendencia refleja el comportamiento de largo plazo y representa el movimiento general de la serie una vez que se han eliminado los otros tres efectos.
La tendencia se define como el comportamiento subyacente de largo plazo en un fenómeno temporal. Normalmente, se realizan ajustes de funciones matemáticas para representar los diferentes tipos de tendencia.
Supongamos que estamos interesados en evaluar el comportamiento del Índice Nacional de Precios al Consumidor (INPC). El INPC se encuentra en el Banco de Información Económica del INEGI1 y su identificador es 628194.
inegi_id <- "628194"
INPC <- inegi_series(inegi_id, token, database = "BIE")
INPC <- INPC[1:367,]
INPC <- ts_invert(INPC)
INPC <- ts(INPC, frequency = 12, start = c(1994,1))
tail(INPC)## [1] 133.681 134.065 134.336 134.087 134.594 136.003
Como primer paso, observamos el fenómeno:
Como se puede apreciar, hay una línea con pocas variaciones. De manera tal que podemos plantear un modelo matemático que incluya la tendencia lineal \(T_t\) y el componente irregular \(I_t\).
\[ Y_t = f(t_t, I_t)\] Si el modelo fuera aditivo, tomaría la siguiente forma
\[Y_t= T_t+I_t\] En particular, el componente de tendencia puede ser modelado como una ecuación lineal con \(\alpha\) como intercepto y \(\beta\) como pendiente
\[T_t= \alpha + \beta t \] En el contexto de las series de tiempo, el componente irregular se identifica con el término de perturbación estócastica, se llama innovación y su nomenclatura es \(a_t\)
\[Y_t =\alpha + \beta t+a_t \]
Para modelar la tendencia, es necesario construir la variable tiempo
\[y_t=\alpha+\beta t+\epsilon_t \] Para modelar la tendencia es necesario construir un modelo econométrico para la función lineal propuesta
\[\hat{y_t}=\hat{\alpha}+\hat{\beta} t \]
##
## Call:
## lm(formula = INPC ~ tiempo)
##
## Coefficients:
## (Intercept) tiempo
## -6856.630 3.448
El resumen estadístico del modelo se imprime con la función summary()
##
## Call:
## lm(formula = INPC ~ tiempo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.868 -3.113 -1.191 2.739 12.314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.857e+03 4.961e+01 -138.2 <2e-16 ***
## tiempo 3.448e+00 2.469e-02 139.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.176 on 365 degrees of freedom
## Multiple R-squared: 0.9816, Adjusted R-squared: 0.9816
## F-statistic: 1.95e+04 on 1 and 365 DF, p-value: < 2.2e-16
A partir del resumen estadístico del modelo de tendencia lineal se pueden almacenar los coeficientes de la regresión.
## [1] -6856.63
## [1] 3.447922
Con base en dichos coeficientes, se realiza un ajuste con los datos observados y los datos estimados por medio de los coeficientes del modelo lineal.
El análisis de series de tiempo tiene dos grandes utilidades, por un lado es factible evaluar el comportamiento histórico del fenómeno en cuestión para analizar el comportamiento a lo largo del tiempo. Por otro lado, con base en la reproducción de los patrones subyacentes, es factible construir, si existe estabilidad paramétrica en los patrones subyacentes, un comportamiento que pemrita reproducir dichos patrones hacia el futuro. Por medio de la librería forecast se puede alcanzar este objetivo.
Construimos la tendencia lineal con los coeficientes y la trasnformamos a serie de tiempo con la temporalidad y la frecuencia de la variable observada y construimos el prónóstico para un horizonte temporal de cinco momentos en el futuro.
Al imprimir el pronóstico se precia el valor estimado del valor futuro, así como sus intervalos con 80% y 90% de confianza.
tendenciaL <- ts(a+b*tiempo, frequency = 12, start = c(1994,1))
pronostico <- forecast(tendenciaL, h=5)
pronostico## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 2024 123.9759 123.9759 123.9759 123.9759 123.9759
## Sep 2024 124.2633 124.2633 124.2633 124.2633 124.2633
## Oct 2024 124.5506 124.5506 124.5506 124.5506 124.5506
## Nov 2024 124.8379 124.8379 124.8379 124.8379 124.8379
## Dec 2024 125.1253 125.1253 125.1253 125.1253 125.1253
Adicionalmente se puede graficar el pronóstico de la tendencia junto con la tendencia histórica
Una manera mucho mas atractiva de mostrar los resultados es por medio de la función autoplot() de la librería forecast y las opciones gráficas de la librería ggplot2.
library(forecast)
library(ggplot2)
inpc <- window(INPC,start=c(1994,1))
fit1 <- tendenciaL
autoplot(inpc) +
autolayer(pronostico, series="Pronóstico")+
autolayer(fit1, series="Tendencia lineal") +
xlab("Año") +
ylab("INPC (índice)") +
ggtitle("Indice Nacional de Precios al Consumidor") +
guides(colour=guide_legend(title="Pronóstico"))Del modelo se desprende que, caeteris paribus, el incremento mensual del Índice Nacional de Precios al Consumidor es de 3.4479223 y cuenta con un intercepto de -6856.6300375
Modelar una tendencia que , en principio, parece ser cuadrática. La variable es la Tasa de condiciones críticas de ocupación (TCCO) (Porcentaje), cuyo identificador es 444608. La tasa de condiciones críticas de ocupación representa el segmento de la población económicamente activa que cuenta con alguna de las siguientes características:
inegi_id <- "444608"
TCCO <- inegi_series(inegi_id, token, database = "BIE")
TCCO <- ts_invert(TCCO)
TCCO <- ts(TCCO, frequency = 12, start = c(2005,1))
ts.plot(TCCO)\[y_t=\beta_0+\beta_1 t+ \beta_3 t^2 +\epsilon_t\]
tiempo <- time(TCCO)
tiempo2 <- tiempo^2
mod1 <- tslm(TCCO~tiempo+tiempo2)
mod1a<- tslm(TCCO~trend+I(trend^2))
b0<-coef(mod1)[[1]]
b1<-coef(mod1)[[2]]
b2<-coef(mod1)[[3]]
tendenciaC <- b0+b1*tiempo+b2*tiempo2
plot(tendenciaC, col="aquamarine2")Tasa de ocupación parcial y desocupación 1 (TOPD1) (Porcentaje), con identificador 444604. La tasa de ocupación parcial y desocupación es el porcentaje de la población económicamente activa que se encuentra desocupada, más la ocupada que trabajó menos de 15 horas en la semana.
inegi_id <- "444604"
TOPD1 <- inegi_series(inegi_id, token, database = "BIE")
TOPD1 <- ts_invert(TOPD1)
TOPD1 <- ts(TOPD1, frequency = 12, start = c(2005,1))
ts.plot(TOPD1)tiempo <- seq(1, length(TOPD1), 1)
tiempo2 <- tiempo^2
tiempo3 <- tiempo^3
mod3 <- lm(TOPD1~tiempo+tiempo2+tiempo3)
mod3a <- tslm(TOPD1~trend++I(trend^2)++I(trend^3))
summary(mod3)##
## Call:
## lm(formula = TOPD1 ~ tiempo + tiempo2 + tiempo3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7320 -0.6310 -0.1026 0.3904 5.4030
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.481e+00 2.519e-01 33.662 < 2e-16 ***
## tiempo 7.367e-02 9.110e-03 8.087 3.30e-14 ***
## tiempo2 -6.169e-04 8.847e-05 -6.973 3.14e-11 ***
## tiempo3 1.363e-06 2.434e-07 5.600 5.98e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9564 on 234 degrees of freedom
## Multiple R-squared: 0.3492, Adjusted R-squared: 0.3409
## F-statistic: 41.85 on 3 and 234 DF, p-value: < 2.2e-16
b0 <- coef(mod3)[[1]]
b1 <- coef(mod3)[[2]]
b2 <- coef(mod3)[[3]]
b3 <- coef(mod3)[[4]]
tendencia <- b0+b1*tiempo+b2*tiempo2+b3*tiempo3
ts.plot(TOPD1, tendencia, col=1:2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 239 9.453163 9.452620 9.453705 9.452333 9.453992
## 240 9.463922 9.462744 9.465100 9.462120 9.465724
## 241 9.474466 9.472527 9.476405 9.471500 9.477432
## 242 9.484800 9.481996 9.487603 9.480512 9.489087
## 243 9.494926 9.491170 9.498682 9.489182 9.500671
Tasa de presión general (TPRG) (Porcentaje), 444605 polinomial
inegi_id <- "444605"
TPRG <- inegi_series(inegi_id, token, database = "BIE")
TPRG <- ts_invert(TPRG)
TPRG <- ts(TPRG, frequency = 12, start = c(2005,1))
plot(TPRG, type="l")tiempo <- seq(1, length(TPRG), 1)
tiempo2 <- tiempo^2
tiempo3 <- tiempo^3
tiempo4 <- tiempo^4
mod4 <- lm(TPRG~tiempo+tiempo2+tiempo3+tiempo4)
mod4a <- tslm(TPRG~trend+I(trend^2)+I(trend^3)+I(trend^4))
summary(mod4)##
## Call:
## lm(formula = TPRG ~ tiempo + tiempo2 + tiempo3 + tiempo4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.69922 -0.72949 -0.03688 0.59556 2.43773
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.725e+00 2.864e-01 23.485 < 2e-16 ***
## tiempo 5.863e-02 1.653e-02 3.548 0.00047 ***
## tiempo2 -5.334e-04 2.804e-04 -1.903 0.05833 .
## tiempo3 1.343e-06 1.761e-06 0.762 0.44659
## tiempo4 -7.905e-10 3.655e-09 -0.216 0.82896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8613 on 233 degrees of freedom
## Multiple R-squared: 0.4921, Adjusted R-squared: 0.4834
## F-statistic: 56.44 on 4 and 233 DF, p-value: < 2.2e-16
b0 <- coef(mod4)[[1]]
b1 <- coef(mod4)[[2]]
b2 <- coef(mod4)[[3]]
b3 <- coef(mod4)[[4]]
b4 <- coef(mod4)[[5]]
tendencia <- b0+b1*tiempo+b2*tiempo2+b3*tiempo3+b4*tiempo4
ts.plot(TPRG, tendencia, col=1:2)Finanzas públicas > Ingresos y egresos brutos por entidad federativa > Nacional 31 entidades federativas (no incluye Ciudad de México) > Ingresos > Impuestos Clave:441566
inegi_id <- "441566"
impuestos <- inegi_series(inegi_id, token, database = "BIE")
impuestos <- ts_invert(impuestos)
impuestos <- ts(impuestos, start = 1989)
ts.plot(impuestos )\[y_t=T_tI_t \] \[y_t=e^{\beta_0+\beta_1t}e^{a_t} \] \[y_t=e^{\beta_0+\beta_1t+a_t} \] \[ln(y_t)= \beta_0+\beta_1t+a_t \]
limpuestos<- log(impuestos)
tiempo <- time(impuestos)
mod5 <- lm(limpuestos~tiempo)
mod5a <- tslm(log(impuestos)~trend)b0 <- coef(mod5)[[1]]
b1 <- coef(mod5)[[2]]
sub_tendencia <- b0+b1*tiempo
tendencia <- exp(sub_tendencia)
ts.plot(impuestos, tendencia, col=1:2)## Rows: 364 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): fecha
## dbl (4): diarios, diarios_cdmx, totales, totales_cdmx
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dias <- seq(1,length(mexicoc$totales),1)
totales <- mexicoc$totales
datos <-data.frame(dias, totales)
dias2 <- seq(1,387,1)
totales2 <- mexicoc$totales[1:387]
datos2 <-data.frame(dias2, totales2)## 7.646194e+13 (6.77e-01): par = (5e+05 100 4)
## 4.346055e+13 (1.01e+00): par = (798186.9 106.3715 22.29686)
## 2.116812e+13 (2.98e+00): par = (1029113 135.3097 77.76)
## 1.444033e+13 (4.24e+00): par = (1812924 245.4666 158.3111)
## 9.841359e+12 (3.32e+00): par = (3149755 333.1296 170.8134)
## 8.787634e+11 (2.82e-01): par = (3653680 320.6704 175.9763)
## 8.210666e+11 (4.68e-02): par = (3248955 302.4053 159.7596)
## 8.193896e+11 (7.27e-03): par = (3339346 306.6954 163.1357)
## 8.193508e+11 (1.16e-03): par = (3336786 306.5205 162.8819)
## 8.193499e+11 (2.11e-04): par = (3339425 306.6546 162.992)
## 8.193499e+11 (3.89e-05): par = (3339163 306.641 162.9785)
## 8.193499e+11 (7.22e-06): par = (3339238 306.6448 162.9818)
sajuste<- summary(ajuste)
a<- sajuste$coefficients[1]
b<- sajuste$coefficients[2]
c<- sajuste$coefficients[3]
c(a,b,c) ## [1] 3339237.8518 306.6448 162.9818
f<-seq(1,1000,1)
contagiados<-a*exp(-exp(-((f-b)/c)))
#format(contagiados, scientific=F)
options(scipen=999)
contagiados <- ceiling(contagiados)
plot(f,contagiados, type = "h",
xlab = "Días", ylab = "Contagios acumulados", col="blue")## dias2 contagiados
## 995 995 3290684
## 996 996 3290979
## 997 997 3291272
## 998 998 3291564
## 999 999 3291853
## 1000 1000 3292141
g <- data.frame(f, contagiados)
dias3 <- d<-seq(1,length(mexicoc),1)
totales3 <- c(contagiados[1:length(mexicoc)])
datos3 <-data.frame(dias3, totales3)library(ggplot2)
ggplot(g, aes(x = f, y = contagiados, colour = "Curva Gompertz") ) +
geom_line()+ theme(legend.position = "bottom")+
geom_point(data = datos, aes(x=dias, y = totales), size = 1, color="orange")+ theme(legend.position = "bottom")+
#geom_line(data=deriv, aes(x=f1, y=derivada, color="brown"))+
xlab("Días desde el pacientel 0") + ylab("Número de contagios confirmados")+
ggtitle("Modelo de predicción de contagios confirmados de Covid-19.")+
labs(subtitle = "Ajuste de una función de crecimiento tipo Gompertz",
caption = "Elaborado por el Laboratorio de Desigualdad Socioeconómica Regional DeSERLab, UAM-I.")+
geom_hline(yintercept=3300000, color="red", show.legend =TRUE, linetype="dashed")+
geom_vline(xintercept=1000, color="red", show.legend=TRUE, linetype="dashed")+
annotate("text", x=250, y=3100000, label="Valor asintótico: 10,965,030")+
annotate("text", x=600, y=2000000, label="2024-02-19")+
annotate("text", x=650, y=2200000, label="Fecha probable del Plateau")Comportamiento que se presenta de manera regular y sistemática
737124 112 Cría y explotación de animales (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades primarias estacionalidad aditiva
inegi_id <- "737124"
IVFEA <- inegi_series(inegi_id, token, database = "BIE")
IVFEA <- ts_invert(IVFEA)
IVFEA <- ts(IVFEA, frequency = 12, start = c(1993,1))
ts.plot(IVFEA)Seasonal Plot:
library(ggplot2)
ggseasonplot(IVFEA, year.labels=TRUE, year.labels.left=TRUE) +
ylab("Índice base 2018=100") +
ggtitle("Seasonal plot: Cría y explotación de animales")+
xlab("Mes")Polar Seasonal Plot:
ggseasonplot(IVFEA, polar=TRUE) +
ylab("Índice base 2018=100") +
ggtitle("Polar seasonal plot: Cría y explotación de animales")+
xlab("Mes")Subseries Plot:
ggsubseriesplot(IVFEA) +
ylab("Índice base 2018=100") +
ggtitle("Subseries seasonal plot: Cría y explotación de animales")+
xlab("Mes")Lag plots
Here the colours indicate the quarter of the variable on the vertical axis. The lines connect points in chronological order. The relationship is strongly positive at lags 4 and 8, reflecting the strong seasonality in the data. The negative relationship seen for lags 2 and 6 occurs because peaks (in Q4) are plotted against troughs (in Q2)
The window() function used here is very useful when extracting a portion of a time series. In this case, we have extracted the data from ausbeer, beginning in 1992.
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## [1,] 1 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0 0
737123 111 Agricultura (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades primarias estacionalidad multiplicativa
inegi_id <- "737123"
IVFAgr <- inegi_series(inegi_id, token, database = "BIE")
IVFAgr <- ts_invert(IVFAgr)
IVFAgr <- ts(IVFAgr, frequency = 12, start = c(1993,1))
ts.plot(IVFAgr)tiempo<- time(IVFAgr)
dummies <- seasonaldummy(IVFAgr)
mod7 <- lm(IVFAgr~tiempo+dummies)
summary(mod7)##
## Call:
## lm(formula = IVFAgr ~ tiempo + dummies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.948 -5.867 -0.866 5.515 33.595
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3186.91639 105.41709 -30.231 < 0.0000000000000002 ***
## tiempo 1.63915 0.05247 31.241 < 0.0000000000000002 ***
## dummiesJan -19.03954 2.36509 -8.050 0.0000000000000116 ***
## dummiesFeb -34.21718 2.36505 -14.468 < 0.0000000000000002 ***
## dummiesMar -42.03273 2.36502 -17.773 < 0.0000000000000002 ***
## dummiesApr -33.09039 2.36500 -13.992 < 0.0000000000000002 ***
## dummiesMay -18.02952 2.36499 -7.624 0.0000000000002122 ***
## dummiesJun -13.54341 2.36499 -5.727 0.0000000213018423 ***
## dummiesJul -24.23038 2.36499 -10.245 < 0.0000000000000002 ***
## dummiesAug -47.34569 2.36500 -20.019 < 0.0000000000000002 ***
## dummiesSep -54.78588 2.36502 -23.165 < 0.0000000000000002 ***
## dummiesOct -43.00850 2.38370 -18.043 < 0.0000000000000002 ***
## dummiesNov -2.80241 2.38369 -1.176 0.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.385 on 368 degrees of freedom
## Multiple R-squared: 0.8558, Adjusted R-squared: 0.8511
## F-statistic: 182 on 12 and 368 DF, p-value: < 0.00000000000000022
##
## Call:
## lm(formula = log(IVFAgr) ~ tiempo + dummies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31782 -0.05770 -0.00164 0.06746 0.35158
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.6805133 1.1257796 -32.582 < 0.0000000000000002 ***
## tiempo 0.0205659 0.0005603 36.704 < 0.0000000000000002 ***
## dummiesJan -0.1968438 0.0252575 -7.793 0.0000000000000674 ***
## dummiesFeb -0.3788068 0.0252571 -14.998 < 0.0000000000000002 ***
## dummiesMar -0.4855111 0.0252568 -19.223 < 0.0000000000000002 ***
## dummiesApr -0.3677764 0.0252565 -14.562 < 0.0000000000000002 ***
## dummiesMay -0.1877426 0.0252564 -7.433 0.0000000000007487 ***
## dummiesJun -0.1308223 0.0252564 -5.180 0.0000003666968376 ***
## dummiesJul -0.2508066 0.0252564 -9.930 < 0.0000000000000002 ***
## dummiesAug -0.5766381 0.0252565 -22.831 < 0.0000000000000002 ***
## dummiesSep -0.7130341 0.0252568 -28.231 < 0.0000000000000002 ***
## dummiesOct -0.5170141 0.0254562 -20.310 < 0.0000000000000002 ***
## dummiesNov -0.0352248 0.0254561 -1.384 0.167
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1002 on 368 degrees of freedom
## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8891
## F-statistic: 254.9 on 12 and 368 DF, p-value: < 0.00000000000000022
lambda= Box-Cox transformation parameter. If lambda=“auto”, then a
transformation is automatically selected using BoxCox.lambda. The
transformation is ignored if NULL. Otherwise, data transformed before
model is estimated.
Moviemiento de mediano plazo
737126 21 Minería (Índice base 2018=100) Indicadores económicos de coyuntura > Indicador global de la actividad económica, base 2018 > Series originales > Índice de volumen físico > Actividades secundarias
inegi_id <- "737126"
IVFMin <- inegi_series(inegi_id, token, database = "BIE")
IVFMin <- ts_invert(IVFMin)
IVFMin <- ts(IVFMin, frequency = 12, start = c(1993,1))
ts.plot(IVFMin)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1993 1 2 3 4 5 6 7 8 9 10 11 12
## Time Series:
## Start = 1993
## End = 2023
## Frequency = 1
## [1] 106.77938 106.63936 104.92355 116.46393 123.60532 126.04125 124.00032
## [8] 127.91148 130.76937 129.10725 133.86573 136.19819 136.93546 135.83764
## [15] 134.25355 128.09959 123.31333 125.66709 126.42954 128.72318 128.72786
## [22] 126.62119 120.72698 116.26661 106.90804 100.00000 94.62836 94.22457
## [29] 96.64103 100.60956 100.74670
778083 -778114 Aguascalientes (Índice de volumen físico 2018 = 100) Indicadores económicos de coyuntura > Indicador trimestral de la actividad económica estatal (ITAEE), base 2018 > Series Originales > Actividades primarias > 111-112 - Agricultura. Cría y explotación de animales > Índice
inegi_id <- "778083"
IVFEAgs <- inegi_series(inegi_id, token, database = "BIE")
IVFEAgs <- ts_invert(IVFEAgs)
IVFEAgs <- ts(IVFEAgs, frequency = 4, start = c(2003,1))
ts.plot(IVFEAgs)## Qtr1 Qtr2 Qtr3 Qtr4
## 2003 -13.799073 -3.184060 11.267782 5.715351
## 2004 -13.799073 -3.184060
## Qtr1 Qtr2 Qtr3 Qtr4
## 2003 NA NA 55.69887 55.74841
## 2004 55.22556 54.89617
## Qtr1 Qtr2 Qtr3 Qtr4
## 2003 NA NA 0.786847 1.302763
## 2004 2.177347 -3.870247
## [1] -13.799073 -3.184060 11.267782 5.715351
dummy <- descomposicion$seasonal
tendencia <- descomposicion$trend
y <- tendencia +dummy
ts.plot(IVFEAgs,y, col=1:2)descomposicion2 <- data.frame(cbind(l=log(IVFEAgs),
t=time(log(IVFEAgs)),
c=cycle(log(IVFEAgs)),
s=factor(descomposicion$seasonal)))
mod8 <- lm(l~t+c+s, data=descomposicion2)
smod8 <- summary(mod8)
yhat <- mod8$fitted.values
ts.plot(IVFEAgs, exp(yhat), col=1:2)SMA= Simple Moving Average EMA= Exponential Moving Average WMA = Weigthed Moving Average DEMA = Differential Exponential Moving Average EVWMA = Elastic Volume Weigthed Moving Average ZLEMA = Zero Lag Exponential Moving Avergae WVMA = Volume Weigthed Moving Average HMA = Diferencia de dos WMA ALMA = Filtro GAussiano
Precio de Apertura Precio de Cierre Precio Máximo Precio Mínimo Volumen
## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast
##
## Attaching package: 'TSA'
## The following object is masked from 'package:readr':
##
## spec
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
## Qtr1 Qtr2 Qtr3 Qtr4
## 2023 104.6562 105.1373 109.6169 109.2868
## 2024 105.3327 105.2148
ts.plot(IVFEAgs, SMA(IVFEAgs, n=10),
SMA(IVFEAgs, n=3),
SMA(IVFEAgs, n=25),
SMA(IVFEAgs, n=40),col=1:5)## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series
Official statistics agencies (such as the US Census Bureau and the Australian Bureau of Statistics) are responsible for a large number of official economic and social time series. These agencies have developed their own decomposition procedures which are used for seasonal adjustment. Most of them use variants of the X-11 method, or the SEATS method, or a combination of the two. These methods are designed specifically to work with quarterly and monthly data, which are the most common series handled by official statistics agencies. They will not handle seasonality of other kinds, such as daily data, or hourly data, or weekly data. We will use the latest implementation of this group of methods known as “X-13ARIMA-SEATS”. For the methods discussed in this section, you will need to have installed the seasonal package in R.
X-11 method The X-11 method originated in the US Census Bureau and was further developed by Statistics Canada. It is based on classical decomposition, but includes many extra steps and features in order to overcome the drawbacks of classical decomposition that were discussed in the previous section. In particular, trend-cycle estimates are available for all observations including the end points, and the seasonal component is allowed to vary slowly over time. X-11 also handles trading day variation, holiday effects and the effects of known predictors. There are methods for both additive and multiplicative decomposition. The process is entirely automatic and tends to be highly robust to outliers and level shifts in the time series. The details of the X-11 method are described in Dagum & Bianconcini (2016).
library(forecast)
lambda <- BoxCox.lambda(IVFEAgs)
IVFEAgs.fit <- ar(BoxCox(IVFEAgs,lambda))
plot(forecast(IVFEAgs,h=20,lambda=lambda))library(ggplot2)
# Using Fourier series for a "ts" object
# K is chosen to minimize the AICc
IVFEAgs.model <- auto.arima(IVFEAgs, xreg=fourier(IVFEAgs,K=2), seasonal=FALSE)
IVFEAgs.fcast <- forecast(IVFEAgs.model, xreg=fourier(IVFEAgs, K=2, h=8))
autoplot(IVFEAgs.fcast) + xlab("Year")## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = IVFEAgs)
##
## Smoothing parameters:
## alpha: 0.05804195
## beta : 0.5067056
## gamma: 0.446458
##
## Coefficients:
## [,1]
## a 106.90612196
## b 0.07312874
## s1 18.50072372
## s2 1.44704963
## s3 -22.30953812
## s4 7.59852503
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2024 Q3 125.47997 118.18408 132.77587 114.32186 136.63809
## 2024 Q4 108.49943 101.17569 115.82317 97.29873 119.70013
## 2025 Q1 84.81597 77.44276 92.18918 73.53962 96.09232
## 2025 Q2 114.79716 107.34712 122.24720 103.40331 126.19101
## 2025 Q3 125.77249 117.14546 134.39952 112.57858 138.96640
## 2025 Q4 108.79194 100.03611 117.54778 95.40104 122.18284
## 2026 Q1 85.10849 76.18705 94.02992 71.46433 98.75264
## 2026 Q2 115.08968 105.96283 124.21653 101.13136 129.04799
## [1] "USD/MXN"
702255 Total industrias manufactureras (Índice base 2018=100) Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índices de remuneraciones medias reales por persona > Índice > Índice de salarios medios reales por obrero
inegi_id <- "702255"
ISMRO <- inegi_series(inegi_id, token, database = "BIE")
ISMRO <- ts_invert(ISMRO)
ISMRO <- ts(ISMRO, frequency = 12, start = c(2007,1))
ts.plot(ISMRO)ISMRO <- window(ISMRO,start=2007)
ISMRO1 <- meanf(ISMRO,h=10)
ISMRO2 <- rwf(ISMRO,h=10)
ISMRO3 <- snaive(ISMRO,h=10)
ISMRO4 <- croston(ISMRO,h=10)
ISMRO5 <- stlf(ISMRO,h=10)
ISMRO6 <- ses(ISMRO,h=10)
ISMRO7 <- holt(ISMRO,h=10)
ISMRO8 <- hw(ISMRO,h=10)
ISMRO9 <- splinef(ISMRO,h=10)
ISMRO10 <- thetaf(ISMRO,h=10)
autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO1, series="Mean", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO2, series="Naïve", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO3, series="Seasonal naïve", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO4, series="Croston´ Method", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO5, series="Loess descomposition", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO6, series="Exponencial Smoothing", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO7, series="Holt", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO8, series="Holt-Winters", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO9, series="cubic SPLine", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))autoplot(window(ISMRO, start=2007)) +
autolayer(ISMRO10, series="Theta Method", PI=FALSE) +
xlab("Mes") + ylab("Índice base 2018=100") +
ggtitle("pronóstico de Índice de salarios medios reales por obrero") +
guides(colour=guide_legend(title="Pronóstico"))702140 - 702160 311 Industria alimentaria (Índice base 2018=100) Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índice de personal ocupado > Índice > Índice de personal ocupado total
inegi_id <- "702140"
IPOAlim <- inegi_series(inegi_id, token, geography = "00", database = "BIE")
IPOAlim <- ts_invert(IPOAlim)
IPOAlim <- ts(IPOAlim, frequency = 12, start = c(2007,1))
ts.plot(IPOAlim)897 Total (Millones de dólares) Indicadores económicos de coyuntura > Balanza comercial de mercancías de México > Series originales > Saldo estacionaria
Indicadores económicos de coyuntura > Tasas de ocupación, desocupación y subocupación (resultados mensuales de la ENOE, 15 años y más) > Series originales > Población total > Población ocupada > Por posición en la ocupación > Trabajadores por cuenta propia Clave:444568
inegi_id <- "444568"
TCP <- inegi_series(inegi_id, token, database = "BIE")
TCP <- ts_invert(TCP)
TCP <- ts(TCP, frequency = 12, start = c(2005,1))
ts.plot(TCP)Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller
## Warning: package 'tseries' was built under R version 4.3.3
## Warning in adf.test(TCP, k = 12): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: TCP
## Dickey-Fuller = -4.1521, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
ARIMA(p,d,q)=ARIMa(1,0,0)
## Series: TCP
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.6286 22.5630
## s.e. 0.0505 0.1118
##
## sigma^2 = 0.4195: log likelihood = -233.58
## AIC=473.16 AICc=473.26 BIC=483.58
## Series: TCP
## ARIMA(0,1,2)(1,0,1)[12]
##
## Coefficients:
## ma1 ma2 sar1 sma1
## -0.3124 -0.4478 0.8255 -0.6774
## s.e. 0.0640 0.0754 0.0873 0.1087
##
## sigma^2 = 0.3874: log likelihood = -223.01
## AIC=456.02 AICc=456.28 BIC=473.36
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01923807 0.6158412 0.4374254 -0.1560695 1.973706 0.6313186
## ACF1
## Training set 0.02425903
Ocupación, empleo y remuneraciones > Tasas de ocupación, desocupación y subocupación (resultados mensuales de la ENOE, 15 años y más) > Urbana, agregado de 32 ciudades > Tasas complementarias > Tasa de ocupación en el sector informal 1 (TOSI1) > Total Clave:444803
Indicadores económicos de coyuntura > Balanza comercial de mercancías de México > Series originales > Saldo > Total Clave:897
library(inegiR)
inegi_id <- "897"
AT <- inegi_series(inegi_id, token, database = "BIE")
AT <- ts_invert(AT)
AT <- ts(AT, frequency = 12, start = c(2005,1))
ts.plot(AT)Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller
## Warning in adf.test(AT, k = 4): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: AT
## Dickey-Fuller = -6.4478, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
ARIMA(p,d,q)=ARIMa(2,1,2)
## Series: AT
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.4406 -0.0130 -0.3134 -0.3352
## s.e. 0.5125 0.0861 0.5099 0.4022
##
## sigma^2 = 2023007: log likelihood = -3513.41
## AIC=7036.82 AICc=7036.97 BIC=7056.84
Indicadores económicos de coyuntura > Índices de precios > Paridades de poder de compra para el Producto Interno Bruto de los países de la OCDE > Producto interno bruto > En moneda nacional a precios corrientes > Estados Unidos Clave:366709
inegi_id <- "366709"
PEA <- inegi_series(inegi_id, token, database = "BIE")
PEA <- ts_invert(PEA)
PEA <- ts(PEA, frequency = 1, start = c(2005,1))
ts.plot(PEA)Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: PEA
## Dickey-Fuller = -0.90883, Lag order = 2, p-value = 0.934
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: diff(PEA)
## Dickey-Fuller = -2.5561, Lag order = 2, p-value = 0.3605
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: diff(diff(PEA))
## Dickey-Fuller = -4.0223, Lag order = 2, p-value = 0.02248
## alternative hypothesis: stationary
ARIMA(p,d,q)=ARIMa(2,1,2)
## Series: PEA
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.5618 0.4350 -0.8380 -0.0308
## s.e. 0.8387 0.8391 0.7497 0.7194
##
## sigma^2 = 245226518374: log likelihood = -362.54
## AIC=735.08 AICc=738.24 BIC=741.18
Indicadores económicos de coyuntura > Encuesta mensual de la industria manufacturera (EMIM). Serie 2018 > Series Originales > Índices de remuneraciones medias reales por hora trabajada > Índice > Índice de remuneraciones medias reales por hora trabajada del personal dependiente de la razón social, considerando utilidades repartidas > 313 Fabricación de insumos textiles y acabado de textiles Clave:702286
inegi_id <- "702286"
AT <- inegi_series(inegi_id, token, database = "BIE")
AT <- ts_invert(AT)
AT <- ts(AT, frequency = 12, start = c(2005,1))
ts.plot(AT)Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: AT
## Dickey-Fuller = -3.4554, Lag order = 12, p-value = 0.04809
## alternative hypothesis: stationary
Indicadores económicos de coyuntura > Producto interno bruto trimestral, base 2018 > Series Originales > Valores a precios de 2018 > Actividades primarias > 11 Agricultura, cría y explotación de animales, aprovechamiento forestal, pesca y caza Clave:735882 Paso 1 Verificar que sea estacionario en niveles
inegi_id <- "735882"
PEA <- inegi_series(inegi_id, token, database = "BIE")
PEA <- ts_invert(PEA)
PEA <- ts(PEA, frequency = 12, start = c(2005,1))
ts.plot(PEA)Verificar que sea estacionaria Contraste Aumentado de Dickey-Fuller
##
## Augmented Dickey-Fuller Test
##
## data: PEA
## Dickey-Fuller = -1.9733, Lag order = 12, p-value = 0.5874
## alternative hypothesis: stationary
## Warning in adf.test(diff(PEA), k = 12): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff(PEA)
## Dickey-Fuller = -4.1943, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
ARIMA(p,d,q)=ARIMa(2,1,2)
## Series: PEA
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.2335 0.7627 -1.7337 0.7396
## s.e. 0.0496 0.0500 0.0559 0.0557
##
## sigma^2 = 1748135974: log likelihood = -2146.49
## AIC=4302.97 AICc=4303.32 BIC=4318.88
## Series: PEA
## ARIMA(1,1,2)(0,1,1)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.9746 0.0928 -0.5984 -0.6109
## s.e. 0.0321 0.0892 0.1029 0.0774
##
## sigma^2 = 682286200: log likelihood = -1924.99
## AIC=3859.98 AICc=3860.35 BIC=3875.54
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 401.279 24849.31 19062.05 -0.07379344 3.19177 0.5409687
## ACF1
## Training set 0.003021428
## date date_shortcut values notes
## 376 1993-06-01 M6 32.30856 NA
## 377 1993-05-01 M5 32.42110 NA
## 378 1993-04-01 M4 31.47853 NA
## 379 1993-03-01 M3 30.87515 NA
## 380 1993-02-01 M2 30.53460 NA
## 381 1993-01-01 M1 31.14023 NA
Paso 1) análisis gráfico paso 2) PRUERBA de Raíces unitarias ADF y
correlogramas Tiene tendencia, estacionalidad
## Warning in adf.test(ener): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ener
## Dickey-Fuller = 0.64895, Lag order = 7, p-value = 0.99
## alternative hypothesis: stationary
Por lo tanto, la variable ener no es estacionaria
Vamos a inspeccionar los correlogramas
paso 3) aplicamos primeras diferencias
## Warning in adf.test(dener): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dener
## Dickey-Fuller = -11.314, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
Modelar parte no estacional
## Series: ener
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.1065 0.1196 0.0408 0.3533
## s.e. 0.2190 0.1026 0.2240 0.0673
##
## sigma^2 = 7.101: log likelihood = -909.88
## AIC=1829.77 AICc=1829.93 BIC=1849.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06092265 2.64727 1.844564 0.1101307 2.631435 0.4375425
## ACF1
## Training set 0.006414525
## Series: ener
## ARIMA(2,1,2)(1,1,0)[12]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1
## -0.2767 -0.0887 0.1794 0.2591 -0.4928
## s.e. 0.4498 0.3082 0.4408 0.2563 0.0456
##
## sigma^2 = 0.0008792: log likelihood = 773.45
## AIC=-1534.9 AICc=-1534.67 BIC=-1511.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02433139 2.040177 1.320867 -0.004510883 1.939994 0.3133182
## ACF1
## Training set -0.09566372
## date date_shortcut values notes
## 15 2008 Year 20442062 NA
## 16 2007 Year 20251027 NA
## 17 2006 Year 19838804 NA
## 18 2005 Year 18929251 NA
## 19 2004 Year 18537508 NA
## 20 2003 Year 17899318 NA
##
## Augmented Dickey-Fuller Test
##
## data: pib
## Dickey-Fuller = -2.4387, Lag order = 2, p-value = 0.4053
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: dpib
## Dickey-Fuller = -3.1207, Lag order = 2, p-value = 0.1454
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: ddpib
## Dickey-Fuller = -3.7084, Lag order = 2, p-value = 0.04225
## alternative hypothesis: stationary
ARIMA(1,0,0) ARIMa(1,0,3) ARIMa(1,1,0) ARIMa(1,2,0)} ARIMA(1,1,3) ARIMA(1,2,3)
library(broom)
m1<- Arima(pib, order=c(1,0,0))
gm1 <- glance(m1)
AIC1 <- gm1$AIC
ts.plot(pib,m1$fitted, col=1:2)m2<- Arima(pib, order=c(1,0,3))
gm2 <- glance(m2)
AIC2 <- gm2$AIC
m3<- Arima(pib, order=c(1,1,0))
gm3 <- glance(m3)
AIC3 <- gm3$AIC
m4<- Arima(pib, order=c(1,1,3))
gm4 <- glance(m4)
AIC4 <- gm4$AIC
m5<- Arima(pib, order=c(1,2,0))
gm5 <- glance(m5)
AIC5 <- gm5$AIC
m6<- Arima(pib, order=c(1,2,3))
gm6 <- glance(m6)
AIC6 <- gm6$AIC
AICS<-data.frame(c("m1", "m2", "m3", "m4", "m5", "m6"),
c(AIC1, AIC2, AIC3, AIC4, AIC5, AIC6))
colnames(AICS) <- c("modelo", "AIC")
AICS## modelo AIC
## 1 m1 609.9071
## 2 m2 615.7220
## 3 m3 575.7077
## 4 m4 579.8303
## 5 m5 553.2977
## 6 m6 550.3682
655561 Ocupación, empleo y remuneraciones > Indicadores de competitividad laboral > Salarios en México por subsector de actividad en la industria manufacturera > Mensual > 336 Fabricación de equipo de transporte
## date date_shortcut values notes
## 76 2018-06-01 M6 2.6 NA
## 77 2018-05-01 M5 2.6 NA
## 78 2018-04-01 M4 2.6 NA
## 79 2018-03-01 M3 2.8 NA
## 80 2018-02-01 M2 2.7 NA
## 81 2018-01-01 M1 2.6 NA
##
## Augmented Dickey-Fuller Test
##
## data: trans
## Dickey-Fuller = -1.384, Lag order = 12, p-value = 0.828
## alternative hypothesis: stationary
##
## Augmented Dickey-Fuller Test
##
## data: dtrans
## Dickey-Fuller = -2.8475, Lag order = 12, p-value = 0.2289
## alternative hypothesis: stationary
## Warning in adf.test(diff(dtrans), k = 12): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff(dtrans)
## Dickey-Fuller = -4.8252, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
SARIMA(4,2,1)(1,0,2) SARIMA(0,2,1)(1,0,2) SARIMA(3,2,1)(1,0,2) SARIMA(2,2,1)(1,0,2)
m1 <- Arima(trans, order = c(0,2,1), seasonal = c(1,1,2))
m2 <- Arima(trans, order = c(3,2,1), seasonal = c(1,1,2))
m3 <- Arima(trans, order = c(2,2,1), seasonal = c(1,1,2))
ts.plot(trans, m1$fitted, m2$fitted,m3$fitted, col=1:4)gm1 <- glance(m1)
AIC1 <- gm1$AIC
gm2 <- glance(m2)
AIC2 <- gm2$AIC
gm3 <- glance(m3)
AIC3 <- gm3$AIC
c(AIC1, AIC2, AIC3)## [1] 135.3360 125.8999 125.1127
We have finished a nice book.
su ubicación específica es Indicadores económicos de coyuntura > Índices de precios > Índice nacional de precios al consumidor > Mensual > Índice↩︎